%--------------------------------------------------------------------------
% SIMULATION
%--------------------------------------------------------------------------


rng('default')

% simulate kappa panel
burn_in         = 1000; % number of initial periods to discart
Nsim            = 100; % number of firms
Tsim            = 120000; % 10,000 years = 120,000 months

[~,states]      = sim_markov_chain_fast(P,Tsim+burn_in,1,1,1:nS);
states          = states';

g_sim           = mu_X(states) + sig_XS(states).*randn(Tsim+burn_in,1) + sig_Xid(states).*randn(Tsim+burn_in,Nsim); % [Tsim,Nsim]

k_def_sim       = k_def(states); % [Tsim,1]

k_sim           = NaN(size(g_sim));
k_sim(1,:)      = k_bar(states(1));

k_QML           = NaN(size(g_sim));
k_QML(1,:)      = k_bar(states(1));

last_issuance   = ones(size(g_sim)); % times initial state of the Markov chain

for t = 1 : Tsim+burn_in-1 % maturity
    
    k_sim(t+1,:) = k_sim(t,:).*exp(g_sim(t,:));
    k_QML(t+1,:) = k_QML(t,:);
    
    default = k_sim(t,:) <= k_def(states(t));
    k_sim(t+1,default) = k_bar(states(t)).*exp(g_sim(t,default));
    k_QML(t+1,default) = k_bar(states(t)); %%% line added 8/4/20 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    issuance = k_sim(t,:) >= k_iss(states(t));
    k_sim(t+1,issuance) = k_bar(states(t)).*exp(g_sim(t,issuance));
    k_QML(t+1,issuance) = k_bar(states(t));
    
    last_issuance(t:end,issuance)  = ones(Tsim+burn_in-t+1,sum(issuance)).*states(t); % most recent issuance state
    last_issuance(t:end,default)  = ones(Tsim+burn_in-t+1,sum(default)).*states(t); % most recent issuance state %%% added 8/4/20 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
end

% eliminate burn_in & and all paths that are already defaulted after burn_in
last_issuance(1:burn_in,:) = [];
k_def_sim(1:burn_in)    = [];
states(1:burn_in)      = [];
k_sim(1:burn_in,:)      = [];
k_QML(1:burn_in,:)      = [];
g_sim(1:burn_in,:)      = [];


default = k_sim <= k_def_sim;


%--------------------------------------------------------------------------
% Table 5
%--------------------------------------------------------------------------


Fit     = griddedInterpolant({1:nS,k},lamDx);
DX_sim  = Fit(last_issuance,k_QML);

Fit     = griddedInterpolant({1:nS,k},lamSx);
SX_sim  = Fit(repmat(states,[1,Nsim]),k_sim);
SX_sim(default) = NaN;

Fit     = griddedInterpolant({1:nS,k},lamS);
S_sim   = Fit(repmat(states,[1,Nsim]),k_sim);
S_sim(default) = NaN;

QMLEV_sim = DX_sim./( DX_sim + SX_sim .* k_sim./k_QML);
QMLEV_sim(default) = NaN;

Fit         = griddedInterpolant({1:nS,k},EP_R);
E_R_sim     = Fit(repmat(states,[1,Nsim]),k_sim);

Fit         = griddedInterpolant({1:nS,k},VP_R');
V_R_sim    = Fit(repmat(states,[1,Nsim]),k_sim);

Fit         = griddedInterpolant({1:nS,k},IV_ATM);
IV_sim      = Fit(repmat(states,[1,Nsim]),k_sim);

Fit        = griddedInterpolant({1:nS,k},IV_Skew);
Skew_sim   = Fit(repmat(states,[1,Nsim]),k_sim);

Fit         = griddedInterpolant({1:nS,k,1:nS},CDS(:,:,:,60));
CDS60_sim   = Fit(repmat(states,[1,Nsim]),k_sim,last_issuance);

ER_sim      = E_R_sim - reshape(Rf(states,1),size(states));

Rcum_sim    = exp(log(S_sim(2:end,:,:))  - log(SX_sim(1:end-1,:,:)) + g_sim(1:end-1,:,:)) - 1;

R_MKT_sim   = trimmean(Rcum_sim',2)'; % expected returns are measured contemporaneously to mktcaps

data = [ QMLEV_sim(:) ER_sim(:) CDS60_sim(:) IV_sim(:) Skew_sim(:) repmat([NaN; R_MKT_sim],[Nsim,1]) V_R_sim(:)];
del = any(isnan(data),2);
data(del,:) = [];

mean_vec    = mean(data(:,1:5));
var_vec     = var(data(:,1:6));
var_vec(2)  = mean(data(:,7)) + var(data(:,2));
table_5_pop = [mean_vec, sqrt(var_vec)]';


%--------------------------------------------------------------------------
% Table 7
%--------------------------------------------------------------------------


Fit          = griddedInterpolant({1:nS,k},PD_q(:,:,60));
PD_Q_60_sim  = Fit(repmat(states,[1,Nsim]),k_sim);

Fit          = griddedInterpolant({1:nS,k},PD_p(:,:,60));
PD_P_60_sim  = Fit(repmat(states,[1,Nsim]),k_sim);

i1 = find(k>k_def(1),1);
i2 = find(k>k_def(2),1);
i3 = find(k>k_def(3),1);
k_small = [k(i1) k(i2) k(i3)]';
states_mat = repmat(states,[1,Nsim]);

Fit  = griddedInterpolant({1:nS,k,1:nS},LGD_q(:,:,:,1));
LGD_Q_1_sim = Fit(states_mat,k_small(states_mat),states_mat);

Fit  = griddedInterpolant({1:nS,k,1:nS},LGD_p(:,:,:,1));
LGD_P_1_sim = Fit(states_mat,k_small(states_mat),states_mat);

data = [PD_P_60_sim(:) PD_Q_60_sim(:) LGD_P_1_sim(:) LGD_Q_1_sim(:)];
data(del,:) = [];

mean_vec    = mean(data);
var_vec     = var(data(:,1:2));
table_7_pop  = [mean_vec(1:2), sqrt(var_vec), mean_vec(3:4)]';


